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Abstract 


In  the  Empirical  Mode  Decomposition  (EMD)  for  the  Hilbert-Huang  Transform  (HHT), 
a  nonlinear  and  nonstationary  signal  is  adaptively  decomposed  by  HHT  into  a  series  of  Intrinsic 
Mode  Functions  (IMFs)  with  the  lowest  one  as  the  trend.  At  each  step  of  the  EMD,  the  low- 
frequency  component  is  mainly  detennined  by  the  average  of  upper  envelope  (consisting  of  local 
maxima)  and  lower  envelopes  (consisting  of  local  minima).  The  high-frequency  component  is 
the  deviation  of  the  signal  relative  to  the  low-frequency  component.  The  fact  that  no  local 
maximum  and  minimum  can  be  determined  at  the  two  end-points  leads  to  detrend  uncertainty, 
and  in  turn  causes  uncertainty  in  HHT.  To  reduce  such  uncertainty,  the  Hermitian  polynomials 
are  used  to  obtain  the  upper  and  lower  envelopes  with  the  first  derivatives  at  the  two  end-points 
(tji,  qR)  as  parameters,  which  are  optimally  detennined  on  the  base  of  minimum  temporal 
variability  of  the  low-frequency  component  in  the  each  step  of  the  decomposition.  This  well- 
posed  mathematical  system  is  called  the  Derivative-optimized  EMD  (DEMD).  With  the  DEMD, 
the  end  effect,  and  detrend  uncertainty  are  drastically  reduced,  and  scales  are  separated  naturally 
without  any  a-priori  subjective  selection  criterion. 
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transform  (HHT),  empirical  mode  decomposition  (EMD);  compact  difference;  Hermitian 
polynomials;  intrinsic  mode  function  (IMF);  end  effect;  detrend  uncertainty 
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1.  Introduction 


Analysis  of  non-stationary  time  series  in  tenns  of  nonlinear  dynamics  has  drawn  attention  in 
many  disciplines.  Traditional  methods,  based  on  linear  and  stationary  assumptions,  are  not 
suitable  to  analyze  nonlinear  and  non-stationary  data.  The  Hilbert-Huang  transform  (HHT)  with 
adaptive  empirical  mode  decomposition  (EMD)  [1]  has  been  developed  to  analyze  and 
nonlinear/non-stationary  data.  Being  adaptive  means  that  the  definition  of  basis  functions  has  to 
be  data-dependent,  not  a-priori  defined  basis  functions  (e.g.,  sinusoidal  functions  in 
linear/stationary  time  series  analysis). 

The  EMD  decomposes  a  nonlinear  and  non-stationary  signal  into  several  intrinsic  mode 
functions  (IMFs)  with  the  lowest  varying  IMF  as  the  trend.  Instantaneous  frequency  can  be 
obtained  by  the  Hilbert-Huang  transform  (HHT),  and  then  the  time-frequency-energy 
distribution  characteristics.  An  IMF  is  a  function  that  must  satisfy  two  conditions  according  to 
the  EMD  algorithm  originally  developed:  (a)  the  difference  between  the  number  of  local  extrema 
and  the  number  of  zero-crossings  must  be  zero  or  one;  (b)  the  running  mean  value  of  the 
envelope  defined  by  the  local  maxima  and  the  envelope  defined  by  the  local  minima  is  zero. 
Average  of  the  upper  and  lower  envelopes  is  treated  as  the  low-frequency  component.  Difference 
of  original  signal  versus  the  low-frequency  component  is  regarded  as  the  high-frequency 
component.  Thus,  accurate  detennination  of  two  envelopes  (i.e.,  one  for  local  maxima  and  the 
other  for  local  minima)  is  crucial  for  the  success  of  the  EMD  in  nonlinear/non-stationary  data 
analysis.  For  a  time  series,  the  interior  extrema  are  easily  identified.  However,  these  extrema  are 
not  enough  to  determine  two  well  behaved  fitting  spline  envelopes  near  the  two  end  points  since 
no  local  maximum/minimum  can  be  identified  there. 
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The  easiest  way  is  to  treat  the  two  end  points  as  “frozen”  points,  i.e.  the  two  end  points 
are  on  both  maximum  and  minimum  envelopes.  Such  a  treatment  makes  the  trend  varying  from 
the  first  end  point  to  the  last  end  point.  Other  practice  is  to  extend  data  points  beyond  the  end 
points  so  as  to  carry  out  the  spline  envelope  fitting  over  and  even  beyond  the  existing  data  range; 
such  as  the  wave  extension  method  [1],  local  straight-line  extension  method  [2],  mirror  or  anti¬ 
mirror  extension  [3],  and  self-similarity  [4],  overlapping  sliding  windows  [5],  and  reject 
segments  close  to  the  endpoints  [2],  While  methods  for  extending  data  vary,  the  essence  of  these 
methods  is  to  predict  data  beyond  the  end  points,  a  dauntingly  difficult  procedure  even  for  linear 
and  stationary  processes.  Since  the  original  signal  only  has  the  extrema  in  the  data  series, 
extending  points  beyond  the  two  end  points  is  not  real.  Therefore  the  data  extension  methods 
will  not  solve  the  problem  no  matter  how  much  efforts  have  been  spent.  Besides,  the  end  error 
may  propagate  from  the  ends  to  the  interior  of  the  data  span  that  would  cause  severe 
deterioration  of  the  IMFs  obtained. 

More  recent  approach  to  deal  with  uncertainties  in  EMD  is  to  use  some  post-processing 
with  allowing  the  sum  of  all  IMFs  different  from  the  original  signal.  For  instance,  the  signal  is 
approximately  represented  by  linear  combination  of  original  IMFs  with  weighting  parameters, 
which  are  determined  using  the  least  square  error  relative  to  the  original  signal.  This  algorithm  is 
called  optimal  EMD  (OEMD)  for  one-dimensional  weights  and  bidirectional  optimal  EMD 
(BOEMD)  for  two-dimensional  weight  matrix  so  as  to  facilitate  approximation  by  window  based 
filtering  [6],  However,  the  OEMD  and  BOEMD  are  limited  by  their  block  based  nature  and  use 
of  adaptive  filters  was  proposed  [7],  The  major  weakness  of  this  type  of  approach  is  the 
inequality  between  the  sum  of  all  IMFs  and  the  original  signal. 


4 


Questions  arise:  Can  the  upper  and  lower  envelopes  be  determined  in  a  systematic  way? 
Can  the  local  maximum  and  minimum  be  objectively  and  optimally  determined  at  the  two  end 
points  without  using  either  extrapolation  or  interpolation  (with  extra  point  beyond  the  end 
point)?  Can  the  sum  of  all  IMFs  always  equal  the  original  signal?  These  problems  will  be  solved 
in  this  study  through  using  the  compact  difference  concepts  [8],  [9],  [10]  with  the  Hennitian 
polynomials.  The  upper  and  lower  envelopes  are  obtained  with  the  first  derivatives  at  the  two 
end-points  {qi,  qfi  as  parameters,  which  are  optimally  detennined  on  the  base  of  minimum 
temporal  variability  of  the  low-frequency  component  in  the  each  step  of  the  decomposition.  This 
method,  called  the  derivative-optimized  EMD  (DEMD),  shows  evident  improvement  in  the 
EMD  analysis.  The  rest  of  the  paper  is  organized  as  follows.  Section-2  introduces  the  classical 
EMD  for  the  Hilbert-Huang  transform  (HHT).  Section-3  describes  the  construction  of  upper  and 
lower  envelopes  using  the  Hennitian  polynomials  with  the  first  derivatives  at  two  end-points,  qi 
(at  / 1 )  and  qR  (at  tfi,  as  tuning  parameters.  Section  4  depicts  optimal  detennination  of  (qu  qfi  on 
the  base  of  minimum  temporal  variability  for  the  low-frequency  component.  Section  5  shows  the 
evaluation.  Section  6  presents  the  conclusions. 

2.  HHT 

The  HHT  has  two  steps.  First,  the  process  of  Empirical  Mode  Decomposition  (EMD) 
reduces  the  time-series  under  analysis  into  components,  known  as  Intrinsic  Mode  Functions 
(IMFs).  Let  a  real  signal  x{t)  be  defined  in  the  time  interval  \t\,  with  two  end  points  x\  =  x(/i), 
and  xN  =  x(tN).  The  EMD  method  is  depicted  as  follows.  First,  the  local  minima  (Xj,  j  =  1, 

2,  . . . ,  J)  and  local  maxima  ( xjmax) ,  k  =  1 ,  2,  . . . ,  K)  of  the  signal  x(/,j  are  identified  with  J  =  K  or 
differing  at  most  by  one.  Second,  interpolation/extrapolation  methods  are  used  to  determine  the 
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upper  and  lower  envelopes  [u(t),  l(t)]  for  t  e.\t\,  fa].  The  mean  of  the  two  envelopes  are 
calculated 

m(t)=[u(t)  + I(t)\/2.  (1) 

The  mean  is  subtracted  from  the  signal,  providing  the  high-frequency  component  (Fig.  1) 

h(t)  =  x(t)  -  m(t),  (2) 

which  is  then  checked  if  it  satisfies  the  above  two  conditions  to  be  an  IMF.  If  yes,  it  is 
considered  as  the  first  IMF  and  denoted 

c(t)  =  h(t).  (3) 

It  is  subtracted  from  the  original  signal  and  the  first  residual, 

KO  =  x(t)  - c(t )  ,  (4) 

is  taken  as  the  new  series  to  continue  the  decomposition.  If  h(t)  is  not  an  IMF,  a  procedure  called 
“sifting  process”  is  applied  as  many  times  as  necessary  to  obtain  an  IMF.  In  the  sifting  process, 
h(t)  is  considered  as  the  new  data,  and  the  same  procedure  applies.  The  IMFs  are  orthogonal,  or 
almost  orthogonal  functions  (mutually  uncorrelated).  This  method  does  not  require  stationarity 
and  linearity  of  the  data  and  is  especially  suitable  for  nonstationary  and  nonlinear  time  series 
analysis.  By  construction,  the  number  of  extrema  decreases  when  going  from  one  residual  to  the 
next;  the  above  algorithm  ends  when  the  residual  has  only  one  extrema,  or  is  constant,  and  in  this 
case  no  more  IMF  can  be  extracted;  the  complete  decomposition  is  then  achieved  in  a  finite 
number  of  steps.  The  signal  x(t)  is  finally  written  as 

x(t)  =  YJcp(t)  +  rl,(t),  (5) 

p= i 

where  cp(t )  is  the  /;-th  IMF  and  />(/)  is  the  residual,  with  no  any  oscillation  (i.e.,  non-existence 
of  both  maximum  and  minimum  envelopes)  and  representing  the  trend. 

The  Hilbert  transform  is  conducted  on  each  IMF  cp(t), 
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zp(t)  =  cp(t)  +  \cp(t). 


(6) 


where  i  =  v— T  and  c At)  is  the  Hilbert  transform  of  cp(t),  represented  by 


where  CP  is  the  Cauchy  principal  value  of  the  integral.  The  complex  variable  zp(t)  in  (6)  can  be 
rewritten  by 

(0  =  cp  {t)  +  i  cp  ( t )  =  ap  (t)  exp[i  9p  (f)],  (8) 

to  obtain  the  instantaneous  amplitude  ap(t)  and  the  instantaneous  phase  function  0  it),  and  the 
instantaneous  frequency  is  calculated  by 

(0p  it)  =  dO pit)  I  dt.  (9) 

Recently,  it  has  been  discovered  that  Hilbert  Transform  has  severe  limitation  on  the  data  for 
instantaneous  frequency  computation.  Different  methods  such  as  direct  quadrature  and 
nonnalized  Hilbert  transform  were  discussed  in  details  [10,  11]. 

The  Hilbert  transfonn  conducted  on  each  IMF,  cPit),  is  called  the  Hilbert-Huang 
transform  (HHT)  [1],  The  key  issue  in  the  EMD  method  for  the  HHT  is  the  accurate 
determination  of  the  upper  and  lower  envelopes  {uPit),  lPit)}.  It  is  noted  the  values  of  {uPit), 
l pit)}  are  only  given  at  the  local  maxima  (x{max>,k  =  1,2,...,  if),  and  local  minima 

iXj,j  =  1,2,...,./ ),  but  unknown  at  the  other  time  instances,  especially  at  the  two  end  points 

t\  and  £n  (Fig.  1).  This  causes  uncertain  in  detennining  cp(t)  with  three  long  recognized 
difficulties:  end-point  effect,  and  detrend  uncertainty. 

3.  Hermitian  Polynomials  for  Upper  and  Lower  Envelopes 
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Let  the  signal  be  discretized  as  (x(f),  i  =  1,  2,  N},  and  let  represent 

occurrence  time  instances  for  either  local  maxima  (t,max, t2max,..., t™x),  or  local  minima 
(hmm,Lmm,...,Cn)  >  as  shown  in  Fig.  2.  Let  em  represent  the  local  maxima  (or  minima)  with  qm  its 
first  derivative  at  the  time  instance rm(m  =  1,2 ,  and  let  Am  =  rm+l -rm,  q  =  (t-rm)/ Am  .  A 
cubic  spline 

Pi  (£ )  =  A  (£)  e,n  +  &  (£)  e„,+l  +  </>3  +  ^4  (£)  ^ 


is  used  to  fill  the  gap  between  two  neighboring  local  maxima  (minima)  between  tam  and  t“H+x  with 
given  values  of  em  and  qm.  This  spline  has  the  following  features 

Pi  (°)  =  *m,  PiO)  =  em+l,  dPi  (°)  / dt  =  qm,  dPi  (!)  / dt  =  <lm+l  •  (H) 

Here,  <f> i(£),  A(A  &(£),  A©  are  defined  by  [5] 

^i(^) = i-3^2 +2<^3,  ^)=?>e-2e, 

(i2) 

Let  the  second  derivative  of  e(t)  be  continuous  at  the  each  local  maxima  (or  minima)  time 
instance,  rm  (m  =  2,3,...,M -1) , 

d2e/  dt2  \{Tm_0)=  d2e /  dt%  J(lm+0).  (13) 


Substitution  of  (10)  into  (13)  leads  to  the  compact  difference  schemes  [8] [9], 


Ac/  ,  H-  2  ( A  ,  +  A  ) Q  A  \Q  ,\  —  3 

nrlm— 1  \  m— 1  m  /  *  m  m—Vl m+l 


^  Ae  ,  Ae  ^ 
A-  —  +  Am_j  ' 


v 


A. 


--I 


A 


Ae  =  e  -  e 

m  m+l  m 


m  J 


(14) 
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where  m  =  2,  3,  M- 1.  Since  the  value  of  e(t)  is  unknown  at  the  two  end-points  (Fig.  2),  the 

condition  (14)  cannot  be  satisfied  at  the  first  and  last  local  maxima  (minima)  time  instances 
xx  and  rM. 

Values  and  derivatives  of  the  local  maximum  (or  minimum)  at  the  two  end-points,  [ e(t\ ), 
e(fv),  e  ’(fv)],  are  unknown.  To  reduce  difficulty,  only  the  first  derivative  at  t\  ( qL )  along 

with  both  value  and  first  derivative  at  x\  (e\,  q\)  are  used  to  fill  the  gap  between  t\  and  x\  with  a 
quadratic  spline, 

Pi^)  =  ¥x^yi+¥2^)qAL+¥^)(lAu  C  =  (t-tx)/AL,  A L=xl-tl,  (15) 

where 

¥,(()  =  1.  ^(0  =  (-l  +  2(-c2)/2,  r,(0  =  (C1-i)/2  (16) 

This  spline  has  the  following  features 

p2(l)  =  ex,  dp2(0) / dt  =  qL,  dp2(l) / dt  =  qx.  (17) 

Continuity  of  the  second  derivative  (13)  at  the  first  local  maximum  (or  minimum)  time  instance 
(A)  from  (t\  -  0)  [using  the  quadratic  spline  (15)]  to  (ti  +  0)  [using  the  cubic  spline  (10)]  leads  to 

— di  (4Ai  +  Aj )  qx  +  2ALq2  —  6AL  —  , 

1  (18) 

A  similar  algebraic  equation  is  obtained  for  the  end-point  tN, 

i  +  (^A^  +  Am1  )  qM  —  A M_xqR  —  6AR  —  -  ,  AR  =  tN  —  xM 

(19) 

where  qR  is  the  first  derivative  at  the  right  end-point.  Fig.  2  shows  the  variation  of  spline  with 
different  values  of  qL  and  qR.  For  simplicity  without  loss  generality,  the  upper  and  lower 
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envelopes  are  assumed  to  have  the  same  first  derivatives  ( qL ,  qR)  at  the  two  end-points  (see  Fig. 


3). 

4.  Optimal  Determination  of  (<//.,  qR) 

Equs.(14),  (18),  and  (19)  represent  two  sets  of  tri-diagonal  linear  equations  for  the  first 
derivatives  at  the  local  maximum  points  of  the  upper  envelope  (quk,  k  =  1.2,3,..., K  )  and  the 

lower  envelope  (<jf,  j  =  1,2,3,...,  J )  as  the  dependent  variables.  Let  the  solutions  for  the  upper 
envelope  be  given  by 

q‘‘  =al+bl‘qL+cukqR,  k  =  1,2,3,...,  K  (20) 

and  the  lower  envelope  be  given  by 

q'j  =  a!j  +  b'jqL  +  c\qR ,  i  =  1,2,3,...,  J.  (21) 

At  each  step,  the  EMD  is  to  decompose  the  signal  into  high  and  low  frequency  components 
with  the  average  of  upper  and  lower  envelopes  using  Eq.(l),  i.e.,  m(t),  as  the  low  frequency 
component  and  the  anomaly  from  the  low  frequency  component  as  the  high  frequency 
component.  Thus,  the  low  frequency  component,  m(t),  should  have  minimum  temporal 
variability.  Usually,  small  absolute  values  of  derivatives  mean  small  temporal  variation.  Since 
the  first  and  second  derivatives  are  already  used  in  obtaining  the  upper  and  lower  envelopes  [see 
Eqs.  (13),  (14)],  minimization  of  integrated  squared  values  of  the  third  derivatives, 

0  r d^u 

S  =  — y  dt  +  — j  dt  — »  mm  ,  (22) 

~  dt  ~  dt 

h  L  -1  H  L  -1 

is  used  to  determine  (qL,  qR).  Here,  the  third  derivative  is  calculated  numerically  by 


k 
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(23) 


VO 

=  6 

V  y 

j  l 

2lj  lj+!  +- 


wr  wr 


yrr  (aJ  +  aJ+i  +  (bj  +  b‘j+l )  ^  +  (CJ  +  c;+i )  c1r  ) 


=  G  XqL,qR) 
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The  end-point  derivatives  (<4,  qR )  are  detennined  by  the  minimization  of  S, 
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As  soon  as  the  end-point  derivatives  (qu  qf)  are  calculated,  the  upper  and  lower  envelopes  are 
detennined.  The  EMD  can  be  effectively  conducted. 

5.  Example 

A  time  series  of  {x,}  (Fig.  4)  consisting  of  a  quadratic  trend  and  three  harmonics, 

3 

*(t, )  =  /o (t, )+Yjfk (*<)>  /o (h )  =  A t] ,  fk (h )  =  A sin (G>ktt +<pk),  k  =  l, 2, 3 

k= 1 

xi=x(tj),  ti  =(i-l)At,  h  =  0,  ^  =  0.9  s,  At  =  0.0018  s,  A  =  501,  (27) 

is  used  to  demonstrate  the  capability  of  DEMD.  The  parameters  in  (27)  are  given  in  Table  1. 
The  left  panels  of  Fig. 5  show  the  trend  and  three  harmonics  of  the  data  represented  by  (27):  fa  = 
=fi(ti),  and  f,  =/3 (tj).  Obviously,  only  /0(f)  represents  the  trend  of  the  signal 
{xi}.  The  trend  of  Jx,}  varies  from f0l  to  f0N, 

/oi  =  0,  f0N=  0.405.  (28) 

The  DEMD  is  conducted  on  the  time  series  {x,-}  (Fig.  4)  to  obtain  three  IMFs  and  a  trend 
(right  panels  of  Fig.  5).  Obviously,  IMFs  well  correspond  to  the  harmonics,  c,(t)  versus  f{t),  with 
high  CCs  (0.978  between  c\  and  f,  0.992  between  c2  and  fo,  0.999  between  c3  and  /3,  1.0 
between  two  trends)  and  low  RRMSEs  (0.00394  between  c\  and  f,  0.0218  between  c2  and^, 
0.00542  between  c3  and  f,  0.0184  between  the  two  trends). 

6.  Conclusions 

(1)  A  major  difficulty  in  HHT  (i.e.,  unknown  local  maximum  and  minimum  at  the  two 
end-points)  has  been  overcome  using  DEMD,  which  is  developed  on  the  base  of  (a)  compact 
difference  scheme  concepts,  and  (b)  minimum  temporal  variability  for  the  low-frequency 
component  (i.e.,  average  of  upper  and  lower  envelopes).  Determination  of  either  upper  or  lower 
envelope  becomes  a  well-posed  mathematical  problem. 
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(2)  The  DEMD  uses  hybrid  Hermitain  polynomials  to  detennine  the  upper  and  lower 
envelopes  with  the  first  derivatives  at  the  two  end-points  (qt,  qR )  as  parameters.  A  set  of  two 
algebraic  equations  for  ( qL ,  qR)  are  derived  on  the  base  of  minimal  integrated  absolute  value  of 
third-order  derivative  (equivalent  to  minimal  temporal  variability).  The  upper  and  lower 
envelopes  are  obtained  after  the  optimal  (qL,  qR)  are  determined. 

(3)  Capability  of  the  DEMD  for  eliminating  end  effect,  and  detrend  uncertainty  is 
demonstrated  using  a  time  series  consisting  of  a  quadratic  trend  and  four  hannonics.  Numerical 
experiment  demonstrates  that  the  new  approach  could  indeed  eliminate  end  effect  and  detrend 
uncertainty  effectively  with  low  RRMSEs.  The  potential  advantage  of  DEMD  over  the 
postprocessing  approaches  is  that  the  sum  of  all  IMFs  is  always  the  same  as  the  original  signal 
using  the  DEMD,  and  is  usually  different  from  the  original  signal  using  the  postprocessing 
approaches. 

(4)  The  three  synthetic  sine  waves  are  used  as  an  example  to  show  the  capability  of  DEMD 
since  exact  components  exist  for  the  error  estimation.  Further  justification  using  real  world 
example  is  needed  although  it  is  difficult  because  exact  components  are  usually  unknown. 
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Table  1.  Values  of  parameters  used  in  Eq.  (27). 
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Amplitude 


Envelopes  and  the  Mean  :  data 


Time  :  second 


1 .  Procedure  of  traditional  EMD. 
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END 


END 


Fig.  3.  Optimal  determination  of  ( qL ,  cjr)  by  minimal  temporal  variability  for  the  upper 
(or  lower)  envelope,  i.e.,  minimal  integrated  square  of  third  derivative  (solid  curves). 
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END 


Original  Signal 


Fig.  4.  Time  series  {jc,-}  represented  by  Eq(27). 
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Components  of  Time  Series 


Time  (s) 

Fig.  5.  Components  of  time  series  {x,}:  a  quadratic  trend  and  four  harmonics  given  by  Eq.(27). 
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Fig.  6.  .  DEMD  on  the  signal  shown  in  Fig.  5:  (a)  IMF-1,  (b)  IMF-2,  (c)  IMF-3,  and  (d)  trend. 
Comparison  between  Fig.  6  and  Fig.  5  shows  the  capability  to  reduce  the  end  effect  and  detrend 
uncertainty. 
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